using EMTSim
using BlockSystems
using ModelingToolkit
using NetworkDynamics
using Graphs
using OrdinaryDiffEq
using DiffEqCallbacks
using SteadyStateDiffEq
using Plots
using Unitful
using CSV
using DataFramesConstants and unit stuff.
ω0 = 2π*50u"rad/s"
Sbase = 300u"MW"
Vbase = 110u"kV" #* sqrt(2/3)
Ibase = Sbase/(Vbase )#* √3) # why the √3 ?
Cbase = Ibase/Vbase
Lbase = Vbase/Ibase
Rbase = (Vbase^2)/Sbase
Rline = 1u"Ω" / Rbase |> u"pu"
Pload = 300u"MW" / Sbase |> u"pu"
Rload = (1u"pu")^2 / Pload |> u"pu" # R=U^2/P
Cline = (2e-6)u"F" / Cbase |> u"s"
Lline = (1/100π)u"H" / Lbase |> u"s"Slack Bus
@variables t u_d(t) u_q(t)
dt = Differential(t)
slackblock = IOBlock([dt(u_d) ~ 0, dt(u_q) ~ 0], [], [u_d, u_q]; name=:slack)IOBlock :slack with 2 eqs
├ inputs: (empty)
├ outputs: u_d(t), u_q(t)
├ istates: (empty)
└ iparams: (empty)create ODE Vertex from this block
slack = ODEVertex(slackblock)R Load
@variables t i_d(t) i_q(t)
@parameters u_d(t) u_q(t) R
loadblock = IOBlock([i_d ~ -1/R * u_d,
i_q ~ -1/R * u_q],
[u_d, u_q], [i_d, i_q],
name=:load)IOBlock :load with 2 eqs
├ inputs: u_d(t), u_q(t)
├ outputs: i_d(t), i_q(t)
├ istates: (empty)
└ iparams: RThe load is used as a current source in a BusBar
busblock = BusBar(loadblock; name=:loadbus)
busblock = set_p(busblock, Dict(:C=>ustrip(u"s", Cline), :ω0=>ustrip(u"rad/s", ω0)))IOBlock :loadbus with 2 eqs
├ inputs: i_d(t), i_q(t)
├ outputs: u_d(t), u_q(t)
├ istates: (empty)
└ iparams: load₊RThe BusBar block can be used to create an ODE Vertex
load = ODEVertex(busblock, [:load₊R])ODE edge for the RL Line
@variables t i_d(t) i_q(t)
@parameters u_d_src(t) u_q_src(t) u_d_dst(t) u_q_dst(t) R L ω
lineblock = IOBlock([dt(i_d) ~ ω * i_q - R/L * i_d + 1/L*(u_d_src - u_d_dst),
dt(i_q) ~ -ω * i_d - R/L * i_q + 1/L*(u_q_src - u_q_dst)],
[u_d_src, u_q_src, u_d_dst, u_q_dst],
[i_d, i_q],
name=:RLLine)
lineblock = set_p(lineblock, Dict(:R=>NoUnits(Rline), :L=>ustrip(u"s", Lline), :ω=>ustrip(u"rad/s", ω0)))IOBlock :RLLine with 2 eqs
├ inputs: u_d_src(t), u_q_src(t), u_d_dst(t), u_q_dst(t)
├ outputs: i_d(t), i_q(t)
├ istates: (empty)
└ iparams: (empty)This block can be used to create an ODEEdge
edge = ODEEdge(lineblock)Set up the network
g = SimpleGraph(2)
add_edge!(g, 1, 2)
nd = network_dynamics([slack, load], edge, g)we use SteadyStateDiffeq to find the steady state
uguess = zeros(length(nd.syms))
uguess[1] = 1.0 # set the d component of the slack to 1 from the beginning
p = ([0, NoUnits(Rload)], nothing)
ssprob = SteadyStateProblem(nd, uguess, p)
u0 = solve(ssprob, DynamicSS(AutoTsit5(Rosenbrock23())))
# SimulationWe simulate a disconnection of the load at t=0.1 s
tspan = (0.0, 0.2)
cb = PresetTimeCallback(0.1, integrator -> integrator.p = ([0, Inf], nothing))
prob = ODEProblem(nd, u0, tspan, p; callback=cb)
sol = solve(prob, Tsit5())
# plot the dq voltage of node 2
# plot(sol, vars=[3,4]; xlims=(0.095, 0.124))retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 1331-element Vector{Float64}:
0.0
0.004596192622990568
0.0048673680168831515
0.005112978790178177
0.005440094588913402
0.005671217534552417
0.0059000150078427115
0.006246626451068669
0.006552133024265822
0.006826198260014092
⋮
0.1989440189891149
0.19909314942295595
0.19923127797379012
0.1993984607587681
0.1995310996834446
0.1996784230803493
0.19983249459047903
0.19996739423890775
0.2
u: 1331-element Vector{Vector{Float64}}:
[1.0, -2.8853683219220016e-22, 0.9758038843594189, -0.024221291204110477, 0.9764177046988668, 0.000507707199198463, -0.9764177046988668, -0.000507707199198463]
[1.0, -2.8853683219220016e-22, 0.9757986924082586, -0.02422368266229332, 0.9764130606766801, 0.0005111425262533142, -0.9764130606766801, -0.0005111425262533142]
[1.0, -2.8853683219220016e-22, 0.9758024273244535, -0.024226587122023014, 0.9764123064894793, 0.0005064159745894732, -0.9764123064894793, -0.0005064159745894732]
[1.0, -2.8853683219220016e-22, 0.9758052654054002, -0.02422284293245541, 0.9764170035367953, 0.0005057270922529847, -0.9764170035367953, -0.0005057270922529847]
[1.0, -2.8853683219220016e-22, 0.9758117316972282, -0.024222076573326895, 0.976420833145115, 0.0005003714301078855, -0.976420833145115, -0.0005003714301078855]
[1.0, -2.8853683219220016e-22, 0.9758045116236415, -0.024220187042279163, 0.9764189884120059, 0.0005076895743561899, -0.9764189884120059, -0.0005076895743561899]
[1.0, -2.8853683219220016e-22, 0.9758037781083057, -0.02422116133412916, 0.9764177679573095, 0.0005078646034726465, -0.9764177679573095, -0.0005078646034726465]
[1.0, -2.8853683219220016e-22, 0.9758029539373411, -0.024221133678456772, 0.9764173908031514, 0.0005086081298502008, -0.9764173908031514, -0.0005086081298502008]
[1.0, -2.8853683219220016e-22, 0.9758020633459226, -0.024222572316079737, 0.9764156824341914, 0.0005086959362062117, -0.9764156824341914, -0.0005086959362062117]
[1.0, -2.8853683219220016e-22, 0.9758036600018949, -0.0242236111649951, 0.976415540326711, 0.0005067753622784793, -0.976415540326711, -0.0005067753622784793]
⋮
[1.0, -2.8853683219220016e-22, 1.000669569462355, -0.0006422832059868946, -1.2912966491950219e-5, 0.02531842394222369, 1.2912966491950219e-5, -0.02531842394222369]
[1.0, -2.8853683219220016e-22, 1.0005880351616798, -0.0006608789007876454, -1.3711426617890149e-5, 0.025384378994485632, 1.3711426617890149e-5, -0.025384378994485632]
[1.0, -2.8853683219220016e-22, 1.0006069611413242, -0.0005980490359880196, 6.192248318034841e-5, 0.02538295439472409, -6.192248318034841e-5, -0.02538295439472409]
[1.0, -2.8853683219220016e-22, 1.000679130244612, -0.0006254352872959244, 9.226484959078049e-6, 0.025317823783546686, -9.226484959078049e-6, -0.025317823783546686]
[1.0, -2.8853683219220016e-22, 1.0006160890452325, -0.0006681106424376321, -3.3780972436427184e-5, 0.02536076536537274, 3.3780972436427184e-5, -0.02536076536537274]
[1.0, -2.8853683219220016e-22, 1.0005851144920301, -0.0006143909445733833, 4.3558069089972444e-5, 0.02539384688969939, -4.3558069089972444e-5, -0.02539384688969939]
[1.0, -2.8853683219220016e-22, 1.0006698769905904, -0.00060301246168973, 4.5377922772656195e-5, 0.025329643442160883, -4.5377922772656195e-5, -0.025329643442160883]
[1.0, -2.8853683219220016e-22, 1.0006510657449301, -0.0006605567589685328, -3.0054199486610306e-5, 0.02533793333598354, 3.0054199486610306e-5, -0.02533793333598354]
[1.0, -2.8853683219220016e-22, 1.000630835548742, -0.0006658911685662713, -3.497205532094123e-5, 0.025352879465425606, 3.497205532094123e-5, -0.025352879465425606]Transformation of the results back to a,b,c frame
a,b,c = Tdqinv(sol.t, sol[3,:], sol[4,:])Comparison with power factory results
Read the Power Facory data and plot for reference
$V_{base}$ is RMS phase-phase voltage.
\[V_{star} = \frac{1}{\sqrt{3}} * V_{base}\\ \hat{V} = \sqrt{\frac{2}{3}} * V_{base}\]
Therefore, i have to multiply the a, b and c results by sqrt(3/2)
df = CSV.read(joinpath(dirname(pathof(EMTSim)), "..", "data","PowerFactory", "Test_EMT.csv"), DataFrame, skipto=3,
header=[:t, :u_1_a, :u_1_b, :u_1_c, :u_2_a, :u_2_b, :u_2_c])
xlims = (0.099,0.124)
plot(sol.t, a.*sqrt(3/2); label="u_2_a", xlims)
plot!(df.t, df.u_2_a; label="PowerFactory A")
plot!(sol.t, b.*sqrt(3/2); label="u_2_b")
plot!(df.t, df.u_2_b; label="PowerFactory B")
plot!(sol.t, c.*sqrt(3/2); label="u_2_c")
plot!(df.t, df.u_2_c; label="PowerFactory C")
... lets have a closer look
xlims!(0.0995,0.105)
This page was generated using Literate.jl.